GLMM Part1

Preparations

Load the necessary libraries

library(car) # for regression diagnostics
library(broom) # for tidy output
library(broom.mixed) ## for tidying mixed effects models
library(ggfortify) # for model diagnostics
library(sjPlot) # for outputs
library(knitr) # for kable
library(effects) # for partial effects plots
library(emmeans) # for estimating marginal means
library(ggeffects) # for partial effects plots
library(MASS) # for glm.nb
library(MuMIn) # for AICc
library(nlme)
library(lme4) # for lmer
library(lmerTest) # for satterthwaite p-values with lmer
library(performance) # for residuals diagnostics
library(see) # for plotting residuals
# library(pbkrtest)  #for kenward-roger p-values with lmer
library(glmmTMB) # for glmmTMB
library(DHARMa) # for residuals and diagnostics
library(tidyverse) # for data wrangling

Scenario

A plant pathologist wanted to examine the effects of two different strengths of tobacco virus on the number of lesions on tobacco leaves. She knew from pilot studies that leaves were inherently very variable in response to the virus. In an attempt to account for this leaf to leaf variability, both treatments were applied to each leaf. Eight individual leaves were divided in half, with half of each leaf inoculated with weak strength virus and the other half inoculated with strong virus. So the leaves were blocks and each treatment was represented once in each block. A completely randomised design would have had 16 leaves, with 8 whole leaves randomly allocated to each treatment.

Tobacco plant

Sampling design

Format of tobacco.csv data files

LEAF	TREAT	NUMBER
1	Strong	35.898
1	Week	25.02
2	Strong	34.118
2	Week	23.167
3	Strong	35.702
3	Week	24.122
...	...	...

LEAF	The blocking factor - Factor B
TREAT	Categorical representation of the strength of the tobacco virus - main factor of interest Factor A
NUMBER	Number of lesions on that part of the tobacco leaf - response variable

Read in the data

tobacco <- read_csv("../public/data/tobacco.csv", trim_ws = TRUE)
glimpse(tobacco)

## Rows: 16
## Columns: 3
## $ LEAF      <chr> "L1", "L1", "L2", "L2", "L3", "L3", "L4", "L4", "L5", "L5", …
## $ TREATMENT <chr> "Strong", "Weak", "Strong", "Weak", "Strong", "Weak", "Stron…
## $ NUMBER    <dbl> 35.89776, 25.01984, 34.11786, 23.16740, 35.70215, 24.12191, …

Exploratory data analysis

The response here represents the average number of lesions across a number of microscope views within a section of leaf. As such, provided the averages are sufficiently far from 0, there is no a priori reason to suspect that the data would not follow a Gaussian distribution.

Model formula: \[ y_i \sim{} \mathcal{N}(\mu_i, \sigma^2)\\ \mu_i =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]

where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of the treatment on the number of lesions. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with leaves.

We will start by explicitly declaring the categorical variable (TREATMENT) as a factor. In addition, random effects (in this case LEAF) should also be declared as factors.

tobacco <- tobacco %>%
  mutate(
    LEAF = factor(LEAF),
    TREATMENT = factor(TREATMENT)
  )

To explore the assumptions of homogeneity of variance and normality, a boxplot of each Treatment level is appropriate.

ggplot(tobacco, aes(y = NUMBER, x = TREATMENT)) +
  geom_boxplot()

Conclusions:

• both normality and homogeneity of variance seem satisfied

It can also be useful to get a sense of the consistency across blocks (LEAF). That is, do all Leaves have a similar baseline level of lesion susceptibility and do they respond similarly to the treatment.

ggplot(tobacco, aes(y = NUMBER, x = as.numeric(LEAF))) +
  geom_line(aes(linetype = TREATMENT))

## If we want to retain the original LEAF labels
ggplot(tobacco, aes(y = NUMBER, x = as.numeric(LEAF))) +
  geom_blank(aes(x = LEAF)) +
  geom_line(aes(linetype = TREATMENT))

Conclusions:

• it is clear that some leaves are more susceptible to lesions (e.g. Leaf 7) than other leaves (e.g. Leaf 4)
• most leaves (other than Leaf 4 and 6) have a similar response to the Treatments - that is most have higher number of lesions from the Strong Treatment than the Weak Treatment.

Given that there are only two levels of Treatment (Strong and Weak), it might be easier to visualise the differences in baselines and effect consistency by plotting as:

ggplot(tobacco, aes(y = NUMBER, x = TREATMENT, group = LEAF)) +
  geom_point() +
  geom_line(aes(x = as.numeric(TREATMENT)))

Conclusions:

• this figure reiterates the points made earlier about the varying baselines and effect consistency.

The above figure also serves as a good way to visualise certain aspects of mixed effects models. When we fit a mixed effects model that includes a random blocking effect (in this case LEAF), we are indicating that we are allowing there to be a different intercept for each block (LEAF). In the current case, the intercept will represent the first Treatment level (Strong). So the random effect is specifying that the intercept can vary from Leaf to Leaf.

We can think of the model as having two tiers (a hierarchy), where the tiers of the hierarchy represent progressively smaller (typically) spatial scales. In the current example, the largest spatial units are the leaves (blocking factor). Within the leaves, there are the two Treatments (Strong and Weak) and within the Treatments are the individual observations.

We tend to represent this hierarchy upside down in the model formula:

\[ y_i \sim{} \mathcal{N}(\mu_i, \sigma^2)\\ \mu_i =\beta_0 + \boldsymbol{\beta} \bf{X_i}\\ \beta_0 = \boldsymbol{\gamma} \bf{Z_i} \]

In addition to allowing there to be a different intercept per leaf, we could allow there to be a different magnitude of effect (difference between Strong and Week Treatment) per leaf. That is considered a random slope. From the figure above, there is some evidence that the effects (slopes) may vary from Leaf to Leaf.

Incorporating a random slope (in addition to a random intercept), may reduce the amount of unexplained variance and thus improve the power of the main effect (Treatment effect).

Fit the model

As with Generalized Linear Models (GLM), Generalized Linear Mixed Models (GLMM) are fitted using maximum likelihood. However, for models containing random effects, maximum likelihood estimates are known to be biased. To correct this bias, rather than use maximum likelihood, we can use Residual Maximum Likelihood (REML). REML maximises the likelihood of the residuals instead of maximising the likelihood of the observed data.

lme (nlme)

## Fit the random intercepts model
tobacco.lme <- lme(NUMBER ~ TREATMENT, random = ~ 1 | LEAF, data = tobacco, method = "REML")
## Fit the random intercepts/slope model
tobacco.lme1 <- lme(NUMBER ~ TREATMENT, random = ~ TREATMENT | LEAF, data = tobacco, method = "REML")

Having now fit the two models (a random intercept model and a random intercept/slope model), we can use AIC (or AICc) to determine which model is ‘best’.

AICc(tobacco.lme, tobacco.lme1)

Conclusions:

• on the basis of AICc, we would conclude that the random intercepts model is the best (it had the lowest AICc). Technically, both AICc values were within 2 units of each other and therefore where not considered significantly different. We would therefore chose models on the basis of complexity (lower complexity, the better). Therefore, we would chose the random intercepts model (if has the fewer used degrees of freedom).

Alternatively, we could use a likelihood ratio test (LRT).

anova(tobacco.lme, tobacco.lme1) %>% print()

##              Model df       AIC       BIC    logLik   Test    L.Ratio p-value
## tobacco.lme      1  4  96.70799  99.26422 -44.35400                          
## tobacco.lme1     2  6 100.67049 104.50484 -44.33525 1 vs 2 0.03749921  0.9814

Conclusions:

• there is no inferential evidence that the two models are different from one another, hence we would go with the simplest model (tobacco.lme - the random intercepts model).

lmer (lme4)

## Fit the random intercepts model
tobacco.lmer <- lmer(NUMBER ~ TREATMENT + (1 | LEAF), data = tobacco, REML = TRUE)
## Fit the random intercepts/slope model
## Note the following could not be run in lmer as there was not enough observations
## to estimate all of the effects and random effects.
## tobacco.lmer1 <- lmer(NUMBER ~ TREATMENT + (TREATMENT|LEAF),  data=tobacco, REML=TRUE)
## tobacco.lmer1 <- lmer(NUMBER ~ TREATMENT + (TREATMENT|LEAF),  data=tobacco, REML=TRUE,
##                       control = lmerControl(
##                           check.nobs.vs.nRE = "ignore",
##                           optimizer = 'optim',
##                           optCtrl = list(method = 'BFGS'))
##                       )

If we had been able to fit both models (the random intercept model and the random intercept/slope model), we could have used AIC (or AICc) to determine which model is ‘best’. As it is, we only have the random intercepts model.

# AICc(tobacco.lme,  tobacco.lme1)

Similarly, we could have used a likelihood ratio test (LRT).

# anova(tobacco.lme,  tobacco.lme1) %>% print

glmmTMB (glmmTMB)

Optimizers include:

• nlminb (non linear optimiser)
• Nelder-Mead - robust but relatively slow
• BFGS - quasi-Newton method (acronym of authors) - gradient based
  • fast, but less robust as requires function to be differentiable
• L-BFGS-B - as above but with different constraints
• SANN - simulated annealing - relatively slow

glmmTMB uses the nlminb (nonlinear minimizer)

## Fit the random intercepts model
tobacco.glmmTMB <- glmmTMB(NUMBER ~ TREATMENT + (1 | LEAF),
  data = tobacco, REML = TRUE
)
## Fit the random intercepts/slope model
## Note the following model did not converge - probably due to insufficient data.
tobacco.glmmTMB1 <- glmmTMB(NUMBER ~ TREATMENT + (TREATMENT | LEAF),
  data = tobacco, REML = TRUE
)
## Try a different optimizer (BFGS)
tobacco.glmmTMB1 <- glmmTMB(NUMBER ~ TREATMENT + (TREATMENT | LEAF),
  data = tobacco, REML = TRUE,
  control = glmmTMBControl(
    optimizer = "optim",
    optArgs = "BFGS"
  )
)
## optArgs = 'Nelder-Mead'))

If we had been able to fit both models (the random intercept model and the random intercept/slope model), we could have used AIC (or AICc) to determine which model is ‘best’. As it is, we only have the random intercepts model.

AICc(tobacco.glmmTMB, tobacco.glmmTMB1)

Similarly, we could have used a likelihood ratio test (LRT).

anova(tobacco.glmmTMB, tobacco.glmmTMB1) %>% print()

## Data: tobacco
## Models:
## tobacco.glmmTMB: NUMBER ~ TREATMENT + (1 | LEAF), zi=~0, disp=~1
## tobacco.glmmTMB1: NUMBER ~ TREATMENT + (TREATMENT | LEAF), zi=~0, disp=~1
##                  Df     AIC     BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## tobacco.glmmTMB   4  96.708  99.798 -44.354   88.708                         
## tobacco.glmmTMB1  6 100.670 105.306 -44.335   88.670 0.0375      2     0.9814

Model validation

lme (nlme)

plot_model

# autoplot(tobacco.lme)
tobacco.lme %>%
  plot_model(type = "diag") %>%
  plot_grid()

Conclusions:

• all the diagnostics look look reasonable

Performance model checking

tobacco.lme %>% performance::check_model()

## Could not compute standard errors from random effects for diagnostic plot.
## Homogeneity of variance could not be computed. Cannot extract residual variance from objects of class 'lme'.

Conclusions:

• all the diagnostics look look reasonable

DHARMa residuals

## unfortunately,  nlme is not supported...
# tobacco.resid = simulateResiduals(tobacco.lme,  plot=TRUE)

lmer (lme4)

plot_model

plot_model() and plot_grid() are of limited value for mixed effects models since:

• the random effects are themselves a list and so do not play nicely with plot_grid()
• the diagnostic plots (such as linearity and homoscadasticity) do not take into account the uncertainty of random effects.

plot_model(tobacco.lmer, type = "diag")[-2] %>% plot_grid()

Conclusions:

• all the diagnostics look look reasonable

Performance model checking

tobacco.lmer %>% performance::check_model()

Conclusions:

• all the diagnostics look look reasonable

DHARMa residuals

tobacco.resid <- tobacco.lmer %>% simulateResiduals(plot = TRUE)

Conclusions:

• all the diagnostics look look reasonable

glmmTMB (glmmTMB)

plot_model

plot_model() and plot_grid() are of limited value for mixed effects models since:

• the random effects are themselves a list and so do not play nicely with plot_grid()
• the diagnostic plots (such as linearity and homoscadasticity) do not take into account the uncertainty of random effects.

plot_model(tobacco.glmmTMB, type = "diag")[-2] %>% plot_grid()

Conclusions:

• all the diagnostics look look reasonable

Performance model checking

tobacco.glmmTMB %>% performance::check_model()

Conclusions:

• all the diagnostics look look reasonable

DHARMa residuals

tobacco.resid <- tobacco.glmmTMB %>% simulateResiduals(plot = TRUE)

Conclusions:

• all the diagnostics look look reasonable

Partial plots

lme (nlme)

plot_model

Unfortunately, plot_model does not support lme models.

# plot_model(tobacco.lme,  type='eff')

allEffects

tobacco.lme %>%
  allEffects() %>%
  plot()

ggpredict

There is a rather obtuse issue with ggpredict. Due to downstream processing (using predict()), models fit with data in the form of a tibble will not work. In order to use ggpredict on a lme model, it is necessary to ensure that the data are data.frame. We can simply update the model..

tobacco.lme1 <- update(tobacco.lme, data = as.data.frame(tobacco))
tobacco.lme1 %>%
  ggpredict(terms = "TREATMENT") %>%
  plot()

ggemmeans

tobacco.lme1 <- update(tobacco.lme, data = as.data.frame(tobacco))
tobacco.lme1 %>%
  ggemmeans(~TREATMENT) %>%
  plot()

lmer (lme4)

plot_model

tobacco.lmer %>% plot_model(type = "eff")

## $TREATMENT

allEffects

tobacco.lmer %>%
  allEffects() %>%
  plot()

ggpredict

tobacco.lmer %>%
  ggpredict(terms = "TREATMENT") %>%
  plot()

ggemmeans

tobacco.lmer %>%
  ggemmeans(~TREATMENT) %>%
  plot()

glmmTMB (glmmTMB)

plot_model

tobacco.glmmTMB %>% plot_model(type = "eff")

## $TREATMENT

allEffects

tobacco.glmmTMB %>%
  allEffects() %>%
  plot()

ggpredict

tobacco.glmmTMB %>%
  ggpredict(terms = "TREATMENT") %>%
  plot()

ggemmeans

tobacco.glmmTMB %>%
  ggemmeans(~TREATMENT) %>%
  plot()

Model investigation / hypothesis testing

lme (nlme)

summary

tobacco.lme %>% summary()

## Linear mixed-effects model fit by REML
##   Data: tobacco 
##        AIC      BIC  logLik
##   96.70799 99.26422 -44.354
## 
## Random effects:
##  Formula: ~1 | LEAF
##         (Intercept) Residual
## StdDev:    3.692213 3.802922
## 
## Fixed effects:  NUMBER ~ TREATMENT 
##                  Value Std.Error DF   t-value p-value
## (Intercept)   34.94013  1.873988  7 18.644797  0.0000
## TREATMENTWeak -7.87938  1.901461  7 -4.143856  0.0043
##  Correlation: 
##               (Intr)
## TREATMENTWeak -0.507
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.61850365 -0.48452438  0.01133404  0.40899997  1.42778194 
## 
## Number of Observations: 16
## Number of Groups: 8

## to get confidence intervals
tobacco.lme %>% intervals()

## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                   lower      est.    upper
## (Intercept)    30.50885 34.940130 39.37141
## TREATMENTWeak -12.37562 -7.879381 -3.38314
## 
##  Random Effects:
##   Level: LEAF 
##                    lower     est.    upper
## sd((Intercept)) 1.580532 3.692213 8.625218
## 
##  Within-group standard error:
##    lower     est.    upper 
## 2.252287 3.802922 6.421125

Conclusions:

• the intercept indicates that the Strong treatment has an average of 34.94 lesions.
• the Weak treatment has 7.88 fewer lesions.
• the variance in intercepts across all Leaves is 13.63
• the scale of variance between Leaves is very similar to the variance within Leaves.

tidy

tobacco.lme %>% tidy(effects = "fixed", conf.int = TRUE)

## including the random effects
tobacco.lme %>% tidy(conf.int = TRUE)

Conclusions:

• the intercept indicates that the Strong treatment has an average of 34.94 lesions.
• the Weak treatment has 7.88 fewer lesions.
• the variance in intercepts across all Leaves is 13.63
• the scale of variance between Leaves is very similar to the variance within Leaves.

tab_model

# warning this is only appropriate for html output
tobacco.lme %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)

	NUMBER
Predictors	Estimates	std. Error	CI	p
(Intercept)	34.94	1.87	30.51 – 39.37	<0.001
TREATMENT [Weak]	-7.88	1.90	-12.38 – -3.38	0.004
Random Effects
σ2	14.46
τ00 LEAF	13.63
N LEAF	8
Observations	16
Marginal R2 / Conditional R2	0.534 / NA
AIC	96.708

\(R^2\)

tobacco.lme %>% r.squaredGLMM()

##            R2m       R2c
## [1,] 0.3707883 0.6761022

## Nakagawa's R2
tobacco.lme %>% performance::r2_nakagawa()

## Random effect variances not available. Returned R2 does not account for random effects.

## # R2 for Mixed Models
## 
##   Conditional R2: NA
##      Marginal R2: 0.534

Conclusions:

• 37% of the total variance in number of lesions is explained by the fixed effects (Treatment) - this is the marginal R2
• 68% of the total variance in number of lesions is explained by the fixed and random effects - this is the conditional R2

lmer (lme4)

summary

tobacco.lmer %>% summary()

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: NUMBER ~ TREATMENT + (1 | LEAF)
##    Data: tobacco
## 
## REML criterion at convergence: 88.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.61850 -0.48453  0.01133  0.40900  1.42778 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  LEAF     (Intercept) 13.63    3.692   
##  Residual             14.46    3.803   
## Number of obs: 16, groups:  LEAF, 8
## 
## Fixed effects:
##               Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)     34.940      1.874 11.332  18.645 7.38e-10 ***
## TREATMENTWeak   -7.879      1.901  7.000  -4.144  0.00433 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## TREATMENTWk -0.507

## to get confidence intervals
tobacco.lmer %>% confint()

##                    2.5 %    97.5 %
## .sig01          0.000000  7.258832
## .sigma          2.333872  6.361451
## (Intercept)    31.216405 38.663855
## TREATMENTWeak -11.829135 -3.929629

Conclusions:

• the intercept indicates that the Strong treatment has an average of 34.94 lesions.
• the Weak treatment has 7.88 fewer lesions.
• the variance in intercepts across all Leaves is 3.8
• the scale of variance between Leaves is very similar to the variance within Leaves.

tidy

tobacco.lmer %>% tidy(effects = "fixed", conf.int = TRUE)

## including the random effects
tobacco.lmer %>% tidy(conf.int = TRUE)

Conclusions:

• the intercept indicates that the Strong treatment has an average of 34.94 lesions.
• the Weak treatment has 7.88 fewer lesions.
• the variance in intercepts across all Leaves is 3.8
• the scale of variance between Leaves is very similar to the variance within Leaves.

tab_model

# warning this is only appropriate for html output
tobacco.lmer %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)

	NUMBER
Predictors	Estimates	std. Error	CI	p
(Intercept)	34.94	1.87	30.86 – 39.02	<0.001
TREATMENT [Weak]	-7.88	1.90	-12.02 – -3.74	0.001
Random Effects
σ2	14.46
τ00 LEAF	13.63
ICC	0.49
N LEAF	8
Observations	16
Marginal R2 / Conditional R2	0.371 / 0.676
AIC	96.708

\(R^2\)

tobacco.lmer %>% r.squaredGLMM()

##            R2m       R2c
## [1,] 0.3707882 0.6761026

## Nakagawa's R2
tobacco.glmmTMB %>% performance::r2_nakagawa()

## # R2 for Mixed Models
## 
##   Conditional R2: 0.676
##      Marginal R2: 0.371

Conclusions:

• 37% of the total variance in number of lesions is explained by the fixed effects (Treatment) - this is the marginal R2
• 68% of the total variance in number of lesions is explained by the fixed and random effects - this is the conditional R2

glmmTMB (glmmTMB)

summary

tobacco.glmmTMB %>% summary()

##  Family: gaussian  ( identity )
## Formula:          NUMBER ~ TREATMENT + (1 | LEAF)
## Data: tobacco
## 
##      AIC      BIC   logLik deviance df.resid 
##     96.7     99.8    -44.4     88.7       14 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  LEAF     (Intercept) 13.63    3.692   
##  Residual             14.46    3.803   
## Number of obs: 16, groups:  LEAF, 8
## 
## Dispersion estimate for gaussian family (sigma^2): 14.5 
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     34.940      1.874  18.645  < 2e-16 ***
## TREATMENTWeak   -7.879      1.901  -4.144 3.42e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tobacco.glmmTMB %>% vcov()

## Conditional model:
##               (Intercept) TREATMENTWeak
## (Intercept)      3.511833     -1.807776
## TREATMENTWeak   -1.807776      3.615552

## the following is not compatible with piping
cov2cor(vcov(tobacco.glmmTMB)$cond)

##               (Intercept) TREATMENTWeak
## (Intercept)     1.0000000    -0.5073298
## TREATMENTWeak  -0.5073298     1.0000000

## to get confidence intervals
tobacco.glmmTMB %>% confint()

##                                    2.5 %    97.5 %  Estimate
## cond.(Intercept)               31.267180 38.613080 34.940130
## cond.TREATMENTWeak            -11.606175 -4.152587 -7.879381
## LEAF.cond.Std.Dev.(Intercept)   1.580539  8.625194  3.692216

tobacco.glmmTMB %>% ranef()

## $LEAF
##    (Intercept)
## L1  -0.3539126
## L2  -1.5406161
## L3  -0.7111775
## L4  -2.5604485
## L5  -1.2399299
## L6   3.5023104
## L7   5.6216063
## L8  -2.7178321

Conclusions:

• the intercept indicates that the Strong treatment has an average of 34.94 lesions.
• the Weak treatment has 7.88 fewer lesions.
• the variance in intercepts across all Leaves is 3.8
• the scale of variance between Leaves is very similar to the variance within Leaves.

tidy

tobacco.glmmTMB %>% tidy(effects = "fixed", conf.int = TRUE)

## including the random effects
tobacco.glmmTMB %>% tidy(conf.int = TRUE)

Conclusions:

• the intercept indicates that the Strong treatment has an average of 34.94 lesions.
• the Weak treatment has 7.88 fewer lesions.
• the variance in intercepts across all Leaves is 3.8
• the scale of variance between Leaves is very similar to the variance within Leaves.

tab_model

# warning this is only appropriate for html output
tobacco.glmmTMB %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)

	NUMBER
Predictors	Estimates	std. Error	CI	p
(Intercept)	34.94	1.87	31.27 – 38.61	<0.001
TREATMENT [Weak]	-7.88	1.90	-11.61 – -4.15	<0.001
Random Effects
σ2	14.46
τ00 LEAF	13.63
ICC	0.49
N LEAF	8
Observations	16
Marginal R2 / Conditional R2	0.371 / 0.676
AIC	96.708

\(R^2\)

tobacco.glmmTMB %>% r.squaredGLMM()

##            R2m       R2c
## [1,] 0.3707882 0.6761025

## Nakagawa's R2
tobacco.glmmTMB %>% performance::r2_nakagawa()

## # R2 for Mixed Models
## 
##   Conditional R2: 0.676
##      Marginal R2: 0.371

Conclusions:

• 37% of the total variance in number of lesions is explained by the fixed effects (Treatment) - this is the marginal R2
• 68% of the total variance in number of lesions is explained by the fixed and random effects - this is the conditional R2

Predictions

Given that there was only two levels of a single predictor, the only follow up analysis might be to explore hypotheses around effects of specific magnitudes. For example, we could investigate whether the strong inoculation treatment increases the number of lesions by over a specific amount. For this example, lets explore the likelihood that the Strong treatment could result in an increase in lesions of over 2 (on average) compared to the Weak treatment.

lme (nlme)

glht

tobacco.lme %>%
  multcomp::glht(linfct = c("TREATMENTWeak=-2")) %>%
  summary()

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lme.formula(fixed = NUMBER ~ TREATMENT, data = tobacco, random = ~1 | 
##     LEAF, method = "REML")
## 
## Linear Hypotheses:
##                     Estimate Std. Error z value Pr(>|z|)   
## TREATMENTWeak == -2   -7.879      1.901  -3.092  0.00199 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Conclusions:

• there is evidence that the Strong treatment does result in a more than 2 unit increase in the number of lesions.

hypothesis

## lme not supported
# brms::hypothesis(tobacco.lme, 'TREATMENTWeak < -2')

lmer (lme4)

glht

tobacco.lmer %>%
  multcomp::glht(linfct = c("TREATMENTWeak=-2")) %>%
  summary()

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lmer(formula = NUMBER ~ TREATMENT + (1 | LEAF), data = tobacco, 
##     REML = TRUE)
## 
## Linear Hypotheses:
##                     Estimate Std. Error z value Pr(>|z|)   
## TREATMENTWeak == -2   -7.879      1.901  -3.092  0.00199 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Conclusions:

• there is evidence that the Strong treatment does result in a more than 2 unit increase in the number of lesions.

hypothesis

tobacco.lmer %>% brms::hypothesis("TREATMENTWeak < -2")

## Hypothesis Tests for class :
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (TREATMENTWeak)-(-2) < 0    -5.91      1.87     -8.9    -2.87        499
##   Post.Prob Star
## 1         1    *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Conclusions:

• there is very strong evidence (Post.Prob = 1) that the Strong treatment does result in a more than 2 unit increase in the number of lesions.

glmmTMB (glmmTMB)

glht

tobacco.glmmTMB %>%
  multcomp::glht(linfct = c("TREATMENTWeak=-2")) %>%
  summary()

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: glmmTMB(formula = NUMBER ~ TREATMENT + (1 | LEAF), data = tobacco, 
##     REML = TRUE, ziformula = ~0, dispformula = ~1)
## 
## Linear Hypotheses:
##                     Estimate Std. Error z value Pr(>|z|)   
## TREATMENTWeak == -2   -7.879      1.901  -3.092  0.00199 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Conclusions:

• there is evidence that the Strong treatment does result in a more than 2 unit increase in the number of lesions.

hypothesis

## Does not appear to work for glmmTMB?
# brms::hypothesis(tobacco.glmmTMB, 'TREATMENTWeak < -2')

Summary figures

lme (nlm)

tobacco.lme %>%
  emmeans(~TREATMENT) %>%
  as.data.frame() %>%
  ggplot() +
  geom_pointrange(aes(y = emmean, x = TREATMENT, ymin = lower.CL, ymax = upper.CL))

lmer (lme4)

tobacco.lmer %>%
  emmeans(~TREATMENT) %>%
  as.data.frame() %>%
  ggplot() +
  geom_pointrange(aes(y = emmean, x = TREATMENT, ymin = lower.CL, ymax = upper.CL))

glmmTMB (glmmTMB)

tobacco.glmmTMB %>%
  emmeans(~TREATMENT) %>%
  as.data.frame() %>%
  ggplot() +
  geom_pointrange(aes(y = emmean, x = TREATMENT, ymin = lower.CL, ymax = upper.CL))

References